import numpy as np

def dirichlet_sample(alphas):
    r = np.random.standard_gamma(alphas)
    r /= r.sum(-1).reshape(-1, 1)
    return r

if __name__ == "__main__":
    alphas0 = np.array([[900,100,2,1,1,1,2,5,1,1,1,1],[100,800,100,2,1,4,1,6,3,4,2,3],[1,100,800,100,3,1,1,1,1,2,4,2],[1,1,2,1,1,1,1,1,2,1,2,1],[800,1,2,2,200,1,1,1,1,1,2,3],[1,1,1,3,1,1,1,1,2,1,1,1],[1,1,800,1,2,3,100,100,2,3,1,1],[1,1,1,2,1,1,1,1,1,1,1,1],[1,1,2,3,800,1,2,1,100,100,3,1],[1,1,2,3,1,1,2,1,100,800,100,1],[1,1,2,1,2,0,800,1,1,100,1,100],[1,1,2,1,1,2,1,800,1,1,100,100]])
    alphas1 = np.array([[900,1,2,1,100,1,2,5,1,1,1,1],[800,200,1,2,1,4,1,6,3,4,2,3],[1,800,100,2,3,1,100,1,1,2,4,2],[1,1,2,1,1,1,1,1,2,1,2,1],[100,1,2,2,800,1,1,1,100,1,2,3],[1,1,2,3,1,1,1,1,2,1,1,1],[100,1,100,1,2,3,800,1,2,3,100,1],[1,1,1,2,1,1,1,1,1,1,1,1],[1,1,2,3,100,1,2,900,1,2,3,1],[1,1,2,3,1,1,2,800,200,1,2,1],[1,1,2,1,2,0,100,1,1,800,100,1],[1,1,2,1,1,2,1,100,1,1,800,100]])
    alphas2 = np.array([[100,100,2,1,800,1,2,5,1,1,1,1],[100,800,100,2,1,4,1,6,3,4,2,3],[1,100,1,100,3,1,800,1,1,2,4,2],[1,1,2,1,1,1,1,1,2,1,2,1],[1,1,2,2,200,1,1,1,800,1,2,3],[1,1,2,3,1,1,1,1,2,1,1,1],[1,1,1,1,2,3,100,100,2,3,800,1],[1,1,1,2,1,1,1,1,1,1,1,1],[1,1,2,3,1,1,2,1,900,100,2,3],[1,1,2,3,1,1,2,1,100,800,100,2],[1,1,2,1,2,1,1,1,1,100,800,100],[1,1,2,1,1,2,1,1,1,1,100,900]])
    alphas3 = np.array([[100,800,2,1,100,1,2,5,1,1,1,1],[1,200,800,2,1,4,1,6,3,4,2,3],[1,1,100,800,3,1,100,1,1,2,4,2],[1,1,2,1,1,1,1,1,2,1,2,1],[100,1,2,2,800,1,1,1,100,1,2,3],[1,1,2,3,1,1,1,1,2,1,1,1],[1,1,100,1,2,3,1,800,2,3,100,1],[1,1,1,2,1,1,1,1,1,1,1,1],[1,1,2,3,100,1,2,1,100,800,2,3],[1,1,2,3,1,1,2,1,1,200,800,2],[1,1,2,1,2,1,100,1,2,3,100,800],[1,1,2,1,1,2,1,100,1,1,1,900]])


    transition_probablity0 = dirichlet_sample(alphas0)
    transition_probablity1 = dirichlet_sample(alphas1)
    transition_probablity2 = dirichlet_sample(alphas2)
    transition_probablity3 = dirichlet_sample(alphas3)
    transitionMatrix= np.dstack((transition_probablity0,transition_probablity1,transition_probablity2,transition_probablity3)) #order is i,j,a


#alphas1 = np.array([[900,1,2,0,100,1,2,5,1,0,0,0],[800,200,1,2,0,4,0,6,3,4,0,3],[1,800,100,2,3,0,100,1,0,2,4,0],[0,0,0,0,0,0,0,0,0,0,0,0],[100,1,2,2,800,1,1,1,100,1,2,3],[0,0,0,0,0,0,0,0,0,0,0,0],[100,1,100,1,2,3,800,1,2,3,100,1],[0,0,0,0,0,0,0,0,0,0,0,0],[1,1,2,3,100,1,2,900,1,2,3,1],[1,1,2,3,1,1,2,800,200,1,2,1],[1,1,2,1,2,0,100,1,1,800,100,1],[1,1,2,1,1,2,1,100,1,1,800,100]])
#alphas2=  np.array([[800,200,1,2,3],[800,2,200,1,3],[800,1,2,2,200],[800,1,4,1,200],[800,4,2,1,200]])

#print(transitionMatrix)
P= transitionMatrix
#print(P[:,:,1])

c = np.array([-0.04, -0.04, -0.04,  +1.0,
              -0.04,   0.0, -0.04,  -1.0,
              -0.04, -0.04, -0.04, -0.04])


r_new = np.tile(c, (12,1))
R=np.dstack((r_new,r_new,r_new,r_new))

print(R[11,:,0])
